System dynamics modeling of lake water management under climate change

Lake Urmia, the twentieth largest lake in the world, is the most valuable aquatic ecosystem in Iran. The lake water level has decreased in recent years due to human activities and climate change. Several studies have highlighted the significant roles of climatic and anthropogenic factors on the shrinkage of the lake. Management policies for water resources harvesting must be adopted to adapt to climate change and avoid the consequent problems stemming from the drought affecting Lake Urmia, and rationing must be applied to the upstream water demands. This study analyzes strategies and evaluates their effectiveness in overcoming the Urmia Lake crisis. Specifically, system dynamics analysis was performed for simulating the water volume of Lake Urmia, and the Hadley Centre coupled model was applied to project surface temperature and precipitation for two future periods: 2021–2050 and 2051–2080. Six management scenarios were considered for decreasing the allocation of agricultural water demand corresponding to two options: (1) one-reservoir option (Bukan reservoir only), and (2) six-reservoir option. The net inflow of Urmia Lake was simulated for the two future periods with the IHACRES model and with artificial neural network models under the six management scenarios. The annual average volumes of Lake Urmia would be 30 × 109 and 12 × 109 m3 over the first and second future periods, respectively, without considering the management scenarios. The lake volumes would rise by about 50% and 75% for the first and second periods, respectively under the management scenarios that involve strict protective measures and elimination of the effect of all dams and their reservoirs. Implementing strict measures would increase the annual average lake volume to 21 × 109 m3 in the second period; yet, this volume would be less than the long-term average and strategic volume. The human water use would be completely eliminated under Scenario 6. Nevertheless, Lake Urmia would experience a considerable loss of storage because of drought.


Scientific Reports
| (2022) 12:5828 | https://doi.org/10.1038/s41598-022-09212-x www.nature.com/scientificreports/ lake have become connected to the land as a result of the lake shrinkage. The high salinity of the lake water has extinguished species such as Artemia 13 . In addition to the severe threat on the wildlife by the lake's drying up salt storms have become more frequent and harmful 14 .
Several studies have been reported addressing multifaceted water issues in Urmia Lake. Delavar and Morid 15 simulated water level fluctuations of Lake Urmia with different methods, such as lake water balance, multiple regression equation (MRE), and artificial neural networks (ANNs). Results showed that the ANN method with concomitant use of accumulative inflow, monthly rainfall, and monthly evaporation had the best accuracy and the least sensitivity in simulating lake's water level variations. Hassanzadeh et al. 16 applied a system dynamics (SD) approach for simulating the effect of surface inflow on the declining Lake Urmia water level. The latter authors concluded that because of drought, over-harvesting of water resources, and dam construction in the basin, the lake water level has declined, and a quarter of the lake area morphed into saline land. Their results showed that the construction of four study reservoirs and surface water over-exploitation with their respective shares of 25% and 65% were the main causes for the lake level reduction. Noury et al. 17 simulated the fluctuations of Lake Urmia water level with support vector machine (SVM) and neural wavelet network (NWN). Results showed that the SVM model performed better than the NWN model. Zarghami and Rahmani 18 examined restoration plans for Lake Urmia with the SD method to increase irrigation efficiency and reduced the acreage of irrigated land. Hosseini-Moghari et al. 19 showed that the water use for anthropogenic purposes in the Lake Urmia basin, being the dominant factor for the Lake Urmia crisis, induced 39-43% of the lake water loss (about 8 km 3 during 2003-2013), up to 90% of the groundwater loss, and 39-45% of the lake-inflow reduction. Bozorg-Haddad et al. 20 demonstrated the accuracy of the SD method in simulating the Urmia Lake water balance and evaluated the impacts of the agricultural water demand supplied by surface water on the Lake Urmia water level.
The severe hydrological drought of Lake Urmia prompted Iran's government to create a ten-year program, named "Urmia Lake Restoration Program (ULRP)" in 2013 with a total cost of five billion US dollars 21 . The ULRP committee has entertained the potential of several alternative management scenarios to control the depletion of Urmia Lake 22 . As such, the following was proposed in the ULRP: reducing and controlling the withdrawal of surface water-groundwater resources in the basin, reducing the water consumption in the agricultural sector by 40%, and transporting water from the Zab and Aras Rivers and from the Caspian Sea 23 . Although numerous studies applied a wide range of methods to show the effects of the factors governing Lake Urmia's level, none of them evaluated the relationship between the allocated water to the agricultural sector from the main reservoirs in the basin (under ULRP management scenarios) and the lake level.
This work's objective is to evaluate management scenarios related to the ULRP to overcome the Lake Urmia crisis in two future periods, namely, 2021-2050 and 2051-2080, under climate change. Since climate change is a serious factor in the present pressing conditions of the Urmia Lake, the results of this research are more practical in comparison to previous studies. The SD method coded in the Vensim software is applied to simulate the current and future water balance of Urmia Lake under climate change. Mitigation strategies proposed in the ULRP are assessed under six management scenarios based on the modeling results for the future periods. This study aims to identify effective restoration measures for the policymakers to resolve the Lake Urmia crisis under climate change.

Methodology
System dynamics method. The SD method applies systemic processing to simulate complex non-linear dynamics and feedback. Systemic processing resorts to various tools to simulate complex system behavior and performance 24 . Systems evolve through states, which change with flows. An example of a state variable is water storage in the study of lakes. The SD method simulates changes in system states driven by flows and various feedbacks 25 .
This work employs the SD method to simulate storage change in Lake Urmia in one historical period  and two future periods (2021-2050 and 2051-2080). The lake's water volume is the state variable, which is governed by inflows (precipitation, surface water inflows, and groundwater inflows) and outflows (evaporation, leakage, and surface water outflows). The lake's mass balance equation is expressed as: where S t+1 , S t , I s , and O s denote the lake's storage at time t + 1, the lake's storage at time t, the inflow rate to the lake at time s (units of volume/time), and the outflow rate from the lake at time s (units of volume/time), respectively.
The SD method employs the Euler and Runge Kutta methods for the solution of differential equations. The software STELLA, Vensim, Powersim, and Dynamo feature SD solvers 26 . This work applies the widely-used Vensim software 27 .
Climate change. The data sets needed for modeling Lake Urmia's storage over the two future periods were generated after simulating the lake's water balance during the historical period. HADCM3, a coupled atmosphere-ocean general circulation model's (AOGCM) climate projections were used to generate precipitation and surface temperature projections over the future periods. The AOGCM data at coarse spatial scales were downscaled to the regional scale suitable for lake storage simulation. The commonly used downscaling methods are statistic and dynamic in nature 28,29 . This works applies the delta-change downscaling method, in which monthly temperature and precipitation differences between the future and historical are calculated by 29 : where ∆T t denotes the difference in long-term average temperatures simulated by HADCM3 for the future ( T GCM,fut,t ) and historical ( T GCM,hist,t ) periods in month t (°C); ∆P t represents the difference in long-term average precipitations simulated by HADCM3 for the future ( P GCM,fut,t ) and historical ( P GCM,hist,t ) periods in month t (mm). Then, ∆T t and ∆P t are applied to project the future downscaled data as follows 29 : where T obs,t, and P obs,t denote respectively the observed temperature (°C) and precipitation (mm) in month t in the baseline period; and T t and P t are the downscaled temperature (°C) and precipitation (mm) in month t of the future period, respectively. Delta-change downscaling is a simple yet efficient option when it comes to spatial downscaling of climate change projections (e.g. [30][31][32]. The gist of this method is to replicate the changing patterns that are projected by the atmospheric ocean general circulation models (AOGCMs) to generate the climate change patterns of hydro-climatic variables on a regional scale. As such, one would simply compute the relative changes in the long-term variations of the variable that is projected by the models within the baseline and future timeframes. These relative changing patterns would be applied to the historical data to project the impact of climate change on a local scale.
Rainfall-runoff modeling. The IHACRES (identification of unit hydrographs and component flows from rainfall, evapotranspiration and streamflow) model is herein applied to simulate runoff from precipitation. Ashofteh et al. 33 implemented the IHACRES model to investigate the effects of climate change on reservoir performance in agricultural water supply. Ashofteh et al. 34 evaluated the probability of flood occurrence in future periods with IHACRES. The IHACRES model includes a non-linear loss module and a linear unit hydrograph module. The non-linear loss module converts the observed rainfall into the effective rainfall, after which the linear unit hydrograph module converts the effective rainfall into the simulated streamflow 35 . Here, precipitation r k in time step k is converted to effective precipitation u k through the non-linear loss module employing a catchment wetness index s k : The effective precipitation is converted to the surface runoff in time step k with the linear unit hydrograph module. The parameters of this model can be set through a thorough grid numeric search and trial-and-error. Perhaps, one of the major advantages of the IHACRES model over other commonly-used rainfall-runoff models is its minimal input data requirement (i.e., air temperature and precipitation) 31,35 .
The other alternative for hydrologic simulation is to use data-driven models. Here, the multilayer perceptron (MLP), a variety of the artificial neural network (ANN) method, was also used to simulate runoff. This model consists of an inlet layer, one or several middle (hidden) layer(s), and an output layer. All of the neurons of a layer are connected to the ones in the next layer, forming a network with complete connections. The primary parameters in modeling the neural network of MLP are: (1) the number of neurons in each layer, (2) the number of layers in the network, and (3) the forcing functions. A regular MLP neural network has three layers 36 . The first and the third layers are respectively the system inputs and outputs. The middle layer consists of neurons that perform calculations on the inputs. Choosing the number of layers in a neural network is made by trial and error 37 . From a hydrological simulation standpoint the main idea behind this model is to create a suitable artificial neural network that is capable of accurately converting a set of hydro-climatic variables such as precipitation and temperature as input data into streamflow values. It should be noted that, like most data-driven models, the process of opting for a proper neural network architecture (i.e., selecting the number of layers, number of neurons, and the forcing function) is, for the most part, a trial-and-error procedure.
One must objectively evaluate the performance of the hydrological models in order to opt for the setting of a suitable parameter. The root mean square error (RMSE), coefficient of determination (R 2 ), and mean absolute error (MAE) are herein employed to assess the performance of the rainfall-runoff model. They are respectively calculated as follows: Performance criteria. Various quantitative measures can be used to assess the performance of water resources systems under different strategies. When it comes to water resources planning and management, perhaps, some of the most common performance criteria are the probability-based performance criteria (PBPC) (i.e., reliability, vulnerability, and resiliency) 31,38 . In this context, reliability represents the probability of successful functioning of a system; resiliency measures the probability of successful functioning following a system failure; lastly, vulnerability is the severity of failure during an operation horizon 39,40 . The basic idea behind a performance evaluation attribute is to provide a quantitative measure to describe and assess the performance of a system. In the context of water resources planning and management, these measures have proven time and again that they can be reliable options to evaluate a set of strategic management options objectively (see, e.g. [40][41][42][43] , and 44 , just to name a few).
Operating policy. Any water resources system requires something called the "rule curve, " which determines how water is allocated in a given situation 45 . A common and effective rule curve when it comes to operation of water resource systems is the standard operation policy (SOP). SOP is a simple, and perhaps best-known realtime operation policy in water resources planning and management 46 . The core principle here is to minimize the water shortage at the current time step with no conservation policy (e.g., hedging rules) in place. The SOP, as a standard rule curve, determines how the operator acts to control a system at any given state of a reservoir 47,48 . This rule curve is established as an attempt to balance various water demands including but not limited to flood control, hydropower, water supply, and recreation 49 . A SOP operating system attempts to release water to meet a water demand at the current time, with no regard to the future. Thus, according to the SOP's principle, the decision-makers, first allocate the available water to meet the demand of the stakeholder with the highest priority. After this first water demand is fully satisfied, the available water can be used for the next demand. Such an allocation process continues until no water is available.

Ethics approval. All authors accept all ethical approvals.
Consent to participate. All authors consent to participate.

Consent to publish.
All authors consent to publish.

Case study
The Lake Urmia basin with an area of 51,862 km 2 is located within 35° 40ʹ and 38° 30 ʹ N (latitude) and within 44° 07 ʹ and 47° 53 ʹ E (longitude), as shown in Fig. 1 (part a). Urmia Lake is the largest lake with the richest aquatic ecosystem in Iran. The lake is located in the northwestern region of Iran between West and East Azarbaijan provinces. The average annual precipitation varies from 250 mm in the lowlands to 1,000 mm in the western regions of the basin. The dominant climate Mediterranean with low precipitation during the summer. The Lake Urmia basin is a typical endorheic (or closed) basin, in which all surface water and groundwater drain into the lake. The Zarinerud and Siminerud Rivers are the two largest rivers in the basin and account for about 52% of the total surface water inflow of Lake Urmia. Figure 1 (part b) depicts a schematic of the river drainage network of the Lake Urmia basin. Figure 2 displays the percentages of runoff contributions from all major tributaries of Lake Urmia.
There are seven large active reservoirs (i.e., Bukan, Mahabad, Alavian, Shahrchai, Hasanlou, Nahand, and Sarough), seven under-construction reservoirs, and 10 under-study reservoirs in the Lake Urmia basin. The Sarough reservoir was not considered in this study because it was recently constructed and there are no operational data prior to 2012. Tables 1 and 2 respectively list the characteristics of the six operational reservoirs in the Lake Urmia basin and their inflows, water demands, and water supplies. The six management scenarios under the climate change condition and their control volumes of the lake's water are summarized in Table 3. The said values are there to show, based on different conditions that are depicted in the said table, how much one should reduce the water demand for the agriculture sector, which by far is the most dominant water-demanding sector in the region. It should be noted that the stated conditions are based on the lake's volume in a given month and its relative status to the long-term average of the lake volume and the lake ecological volume in the said month. The general idea being that the more water stresses are imposed on the lake, the more pronounced restriction must be applied to the demand to control and potentially restore the water of the Urmia Lake. The water allocated to the agricultural sector from either the Bukan reservoir (Option 1) or all six reservoirs (Option 2) were determined under each scenario based on the simulation of lake storage. This study has chosen the SOP as the operating policy. Based on the Urmia Lake Restoration Program (ULRP) recommendations these priorities, from the highest to the lowest one, are drinking, ecological, agricultural, and industrial demands.    Table 3. Six management scenarios and their lake water control volumes (adopted from 22,23 ). V t = lake volume in month t (million m 3 ); V ave = long-term average of the lake volume in month t (million m 3 ); and V eco = lake ecological volume in month t (million m 3 ) (i.e., the minimum lake water volume for the survival of living organisms in the Lake Urmia ecosystem).

Management scenarios
Lake water volume condition

Results and discussion
Simulation model for historical data. The SD method performance in simulating the lake volume was evaluated by comparing simulated values with the observed data. Figure 3 shows the comparison of the observed and simulated water volumes of Urmia Lake in 1957-2005. The values of R 2 and RMSE were 0.952 and 18.2 × 10 6 m 3 , respectively. Figure 4 depicts a comparison between the simulated and observed lake water volumes, indicating that most data points are scattered near the 1:1 (45°) line, although the maximum simulated lake volume was slightly smaller than the observed one. It is concluded that the overall performance of the SD method was acceptable based on the calculated results. Rainfall-Runoff modeling. The number of neurons between 1 and 20 was investigated in the ANN model.
Here, precipitation and temperature were considered as the input variables. 12 neurons were chosen as the optimal number by trial and error. Sensitivity analysis was done on the activation function for the input and output layers. The results showed that the sigmoid function had a suitable performance as a non-linear activator   Results of the performance measures of the two models for calibration and validation are listed in Table 4. The ANN model performed better than IHACRES in terms of R 2 , MAE, and RMSE. In calibration the ANN model had a better performance with an increase of 13.5% in R 2 , and 11% and 33% decreases in RMSE and MAE, respectively, in comparison with the IHACRES model. In validation the ANN model had a 3% increase in R 2 , and 5% and 18% reductions in MAE and RMSE, respectively, in comparison with the IHACRES model. Therefore, the ANN model was implemented to simulate the inflow to the Bukan reservoir in the two future periods. Figure 6 shows the inflow to the Bukan reservoir would increase by 1.07 × 10 6 m 3 and 176.33 × 10 6 m 3 in the first and second future periods (compared to the historical data), respectively. According to Fig. 6 in the first future period the increase of inflow to the Bukan reservoir is evident from January through March, October, and December. However, the inflow would be reduced from April through July. The inflow changes would be more significant in the second future period than the first one. The largest inflow in the second future period would be in March, while the peak inflow in the historical period would be in April. It can be concluded that the inflow to the Bukan reservoir would change in the future.
Reservoir water supplies in the future periods. The water supplies from the Bukan reservoir were calculated for the two future periods. Table 5 lists the monthly potable and industrial water demand, agricul-  Table 5. Monthly potable, industrial, and agricultural water demands from the Bukan reservoir in the two future periods (10 6 m 3 ) (adopted from 22,23  www.nature.com/scientificreports/ tural water demand, and Fig. 7 depicts the environmental water demand from the Bukan reservoir in the future periods 22,23 . Evapotranspiration in the agricultural sector and evaporation from the lake would increase due to the increased temperature in both future periods, and the potable and industrial water demand would increase due to the population growth. As expected, the agricultural, environmental, and potable and industrial water demands would increase. Also, the water inflow to the lake under climate change would increase, so the share of the lake from it would increase. However, the evaporation of the lake could be greater than the increased environmental water demand because of the temperature increment. The increase in the potable and industrial water demand would be negligible in comparison with the rise in the water demands of the agricultural and environment sectors. The potable and industrial demands in various months would be nearly constant (about 17 × 10 6 m 3 with a variation of 1 × 10 6 m 3 ). The agricultural water demand would zero in the four cold months of the year when the river inflow was much higher. Therefore, the Bukan reservoir stores water during the cold season for regulating water for agricultural use during the growth season. Among all water demands, the agriculture sector would have the highest water demand from June through August, and its average monthly value equals 500 × 10 6 m 3 . Figure 7 shows considerable variations in the monthly lake's environmental water demands, which are similar to the variations in the inflows of the lake. The environmental water demand for the lake in the second future period (2051-2080) would be larger than the first future period's (2021-2050).
Simulation of Lake Urmia water stage in the future periods. The variables for water allocation of a reservoir in this study include: reservoir inflow, flow downstream the reservoir, and water demands of the agricultural, environmental, potable, and industrial sectors. These variables for the Bukan reservoir were compared with those for the five other operational reservoirs in the Lake Urmia basin. Specifically, the ratios of these variables for the five reservoirs to the corresponding variables for the Bukan reservoir in the historical period were calculated (see Table 6). The ratios presented in Table 6 were used to calculate the inflow and outflows of the five other reservoirs in the future periods. The reservoir inflow and its downstream flow in the future periods were calculated only for the Bukan reservoir.
The simulated water stages of Lake Urmia in the first and second future periods are respectively shown in Figs. 8 and 9. The lake water level in the period 2021-2050 would be lower than its long-term average (1,275.86 m), but also lower than the ecological water level (i.e., the minimum water level for the survival of living organisms in the Lake Urmia ecosystem is 1,274.10 m). Remedial strategies have been proposed to prevent such a crisis in the basin in the Urmia Lake Restoration Program (ULRP). Specifically, management must be implemented to reduce the agricultural water use. Reducing agricultural water allocation is examined in the following section. In the second future period the lake water level would be slightly higher than the corresponding historical average water level determined by the mass balance analysis. Table 6. Ratios of five variables of five reservoirs to the corresponding values for the Bukan reservoir in the historical period (partially adopted from 22,23 ). IN Inflow; AgD Agricultural water demand; PID Potable and industrial water demand; ED Environmental water demand; DF Downstream flow. The performance criteria for the simulated volumes of Lake Urmia in the period of 2021-2050 corresponding to Option 2 are better than those for Option 1. This is an indication Option 2, in which all six reservoirs would be operational, is preferable so long as meeting agricultural demand is the main objective.
The time-based reliability of the lake water volume for Option 2, with a performance threshold of 70% in the first future period, increases from 0 in management Scenario 1 to 2% in management Scenario 6. This 2% increase in reliability is achieved when the agricultural water demand decreases by 100%, while the value of the time-based reliability of the lake water volume is constantly 0 for Option 1. The lake volumetric reliability increases by an average of 10% in Option 1 for the second future period under all six management scenarios if the performance threshold is reduced from 70 to 50%, while the lake time-based reliability increases by an average of 30% in Option 2.
The simulated lake water levels in 2021-2050 would be lower than those in the 2051-2080 period, which means the values of the four performance criteria for the simulated lake water volumes in the second future period would be better than those in the first future period when the same policy is implemented in both future periods. For instance, under Scenario 1 and Option 1 the vulnerability values decrease from 24%, 43%, 56%, and 60% in 2021-2050 (Fig. 10) to 1%, 5%, 16%, and 22% (Fig. 10) in 2051-2080 with respect to the performance thresholds of 50%, 70%, 90%, and 100%, respectively. Comparison of the management scenarios revealed that there was no significant difference in the performances corresponding to management Scenarios 1 through 4, while all performance criteria were improved in management Scenarios 5 and 6. Our results suggested that implementing stricter policies would help improve the conditions of the lake and help its ecosystem.
The Zarinerud sub-basin and Bukan reservoir play central roles on the Lake Urmia hydrology. Future lake studies should focus on the sub-basins and reservoirs that have govern the lake's hydrology 8,50 . Moghadasi et al. 51 showed that agricultural water use should be reduced by 25-35% in East Azerbaijan, and by 15-25% in West Azarbaijan to achieve optimal water allocations in the Lake Urmia basin.

Conclusions
The temperature and precipitation data were extracted for the Lake Urmia basin from the HADCM3 projections under emission Scenario A2. These projections were downscaled with the delta-change method for future periods 2021-2050 and 2051-2080. The inflows to the Bukan reservoir and the downstream flows were simulated with the ANN model. The inflows of five other reservoirs and their downstream flows were calculated using the Bukan reservoir data and the ratios of the values of five reservoirs' variables to the corresponding values for Bukan reservoir located in the Zarrineh sub-basin. The SD method was applied to simulate the water balance of the Lake Urmia basin in the future periods. Several management scenarios were developed to improve the condition of Urmia Lake. The following conclusions were obtained from this study: • The temperature in the Lake Urmia basin is expected to rise in both future periods. The precipitation in region upstream of the Bukan reservoir is expected to decrease by 32% in the first future period (2021-2050) and increase by 0.4% in the second future period (2051-2080). The precipitation is expected to increase by 0.3% in the first period and decrease by 2.38% in the second period in the regions downstream of the Bukan reservoir. • The lake water level in the first future period would be lower than its long-term average and the ecological water level. • The values of the performance criteria of the predicted lake water volumes in the first future period would be lower than those in the second future period under all management scenarios. But the rate of increase in the performance criteria by decreasing the agricultural water allocation (changing from management Scenario 1 to management Scenario 6) in the first future period would be larger than that of the second future period. • Lake Urmia would experience a considerable loss of water because of drought even if the agricultural water allocation is eliminated under some management scenarios.
More in-depth studies could help shed light on the intricacies and challenges of preserving Lake Urmia.

Availability of data and materials
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The codes that support the findings of this study are available from the corresponding author upon reasonable request.